#####################################
## "Backyard politics in Foreign Aid"
## William Christiansen
## Tobias Heinrich
## Timothy Peterson
#####################################

## File that estimates all models
#################################


## Load data, set basics
########################
load("output/data_d1.Rdata")
load("output/data_d2.Rdata")

form_main <- "+ T1 + T2"
form_covs <- "+ Gender + Age + Costs + Education_hi + Education_lo + Ideology + HousingMed + HousingPctIncrease + Transnational + GreatestCountry"
tmp <- model.matrix(object=formula(paste0("Y1 ~ 1", form_main, form_covs)),
                    data=d1)


#############
## Y | T1, T2
#############
out_YT <- array(data=0, dim=c(n_bs, 2, 2, ncol(tmp)),
                dimnames=list(Iterations=NULL, DV=c("Rating", "Thermometer"), 
                              Scenario=c("Increase", "Cut"),
                              Coeffcient=colnames(tmp)))
for(j in 1:2)  ## Loop through Cut/ increase
{
  this_d1 <- subset(d1, Change == ifelse(j == 1, "Increase", "Cut"))
  this_d2 <- subset(d2, Change == ifelse(j == 1, "Increase", "Cut"))
  n_d1 <- nrow(this_d1)
  n_d2 <- nrow(this_d2)
  
  ## DV: Project rating
  for(iter in 1:n_bs)
  {
    ind <- sample(1:n_d1, size=n_d1, replace=TRUE)
    mod <- lm(formula=as.formula(paste0("Y1 ~ 1", form_main)), data=this_d1[ind,])
    out_YT[iter, "Rating", j, 1:length(coef(mod))] <- coef(mod)
  }
  
  ## DV: Feeling thermometer
  mod <- ols(formula=as.formula(paste0("Y2 ~ 1", form_main)), data=this_d2, model=TRUE,
             x=TRUE, y=TRUE)
  mod <- bootcov(fit=mod, cluster=this_d2$ID, B=n_bs, coef.reps=TRUE)   
  out_YT[, "Thermometer", j, 1:ncol(mod$boot.Coef)] <- mod$boot.Coef
}



##################
## M_i | T1, T2, X
##################
out_MT <- array(data=0, dim=c(n_bs, 6, 2, 2, ncol(tmp)),
                dimnames=list(Iterations=NULL, M=paste0("M", 1:6), 
                              Scenario=c("Increase", "Cut"),
                              CovSet=c("Basic", "Expanded"),
                              Coeffcient=colnames(tmp)))
for(i in 1:2)  ## Loop through covariate spec
{
  for(j in 1:2)  ## Loop through Cut/ increase
  {
    this_d1 <- subset(d1, Change == ifelse(j == 1,  "Increase", "Cut"))
    n_d1 <- nrow(this_d1)
    
    for(m in 1:6)
    {
      this_form <- as.formula(paste0("M", m, " ~ 1 + ", form_main, ifelse(i == 1, "", form_covs)))
      for(iter in 1:n_bs)
      {
        ind <- sample(1:n_d1, size=n_d1, replace=TRUE)
        mod <- glm(formula=this_form, data=this_d1[ind,], family=binomial(link='probit'))
        out_MT[iter, m, j, i, 1:length(coef(mod))] <- coef(mod)
      }
    }
  }
}



#############################
## Y | M1, ..., M6, T1, T2, X
#############################
out_YM <- array(data=0, dim=c(n_bs, 2, 2, ncol(tmp)+6),
                dimnames=list(Iterations=NULL, DV=c("Rating", "Thermometer"), 
                              Scenario=c("Increase", "Cut"),
                              Coeffcient=c(colnames(tmp)[1], paste0("M", 1:6),
                                           colnames(tmp)[-1])))
this_form <- paste0(paste(paste0("M", 1:6), collapse=" + "), form_main, form_covs)
for(j in 1:2)  ## Loop through Cut/ increase
{
  this_d1 <- subset(d1, Change == ifelse(j == 1, "Increase", "Cut"))
  this_d2 <- subset(d2, Change == ifelse(j == 1, "Increase", "Cut"))
  n_d1 <- nrow(this_d1)
  n_d2 <- nrow(this_d2)
  
  ## DV: Package rating
  for(iter in 1:n_bs)
  {
    ind <- sample(1:n_d1, size=n_d1, replace=TRUE)
    mod <- lm(formula=as.formula(paste0("Y1 ~ 1 + ", this_form)), data=this_d1[ind,])
    out_YM[iter, "Rating", j, 1:length(coef(mod))] <- coef(mod)
  }
  
  ## DV: Feeling thermometer
  mod <- ols(formula=as.formula(paste0("Y2 ~ 1 + ", this_form)), data=this_d2, model=TRUE,
             x=TRUE, y=TRUE)
  mod <- bootcov(fit=mod, cluster=this_d2$ID, B=n_bs, coef.reps=TRUE)   
  out_YM[, "Thermometer", j, 1:ncol(mod$boot.Coef)] <- mod$boot.Coef
}



###############################
## Do causal mediation analysis
###############################
out_med <- array(NA, dim=c(n_bs, 2, 2, 2, 8),
                 dimnames=list(Iterations=NULL, DV=c("Rating", "Thermometer"), 
                               Treatment=c("T1", "T2"), 
                               Scenario=c("Increase", "Cut"),
                               Effect=c("Total effect", "Direct effect", 
                                        "Indirect effect (M1)", "Indirect effect (M2)",
                                        "Indirect effect (M3)", "Indirect effect (M4)",
                                        "Indirect effect (M5)", "Indirect effect (M6)")))

for(i in 1:2)  ## Loop through DV
{
  for(j in 1:2)  ## Loop through increase/cut
  {
    for(n in 1:2)  ## Loop through treatments
    {
      ## Total effect
      ###############
      TE <- out_YT[, ifelse(i == 1, "Rating", "Thermometer"), j, paste0("T", n)]
      
      ## Indirect effects for M1, ..., M6
      ###################################
      this_d1 <- subset(d1, Change == ifelse(j == 1, "Increase", "Cut"))
      IE <- matrix(0, n_bs, 6)
      for(m in 1:6)
      {
        m_data <- model.matrix(object=as.formula(paste0("M", m, " ~ 1 + ", form_main, form_covs)),
                               data=this_d1)
        m_data_0 <- m_data_1 <- m_data[m_data[, paste0("T", n)] == 1, ]
        m_data_0[, paste0("T", n)] <- 0
        m_pred_1 <- pnorm(m_data_1 %*% t(out_MT[, m, j, 2, ]))
        m_pred_0 <- pnorm(m_data_0 %*% t(out_MT[, m, j, 2, ]))
        m_diff_10 <- m_pred_1 - m_pred_0
        n_obs <- nrow(m_data_0)
        
        for(l in 1:n_obs)
        {
          IE[, m] <- IE[, m] + out_YM[, ifelse(i == 1, "Rating", "Thermometer"), j, paste0("M", m)] * m_diff_10[l, ]
        }
        IE[, m] <- IE[, m] / n_obs
      }
      
      ## Direct effect
      ################
      DE <- TE - rowSums(IE)
      
      out_med[, i, n, j, ] <- cbind(TE, DE, IE)
    }
  }
}

save(out_med, file="output/est_out_med.Rdata")
save(out_YT, file="output/est_out_YT.Rdata")
save(out_MT, file="output/est_out_MT.Rdata")
save(out_YM, file="output/est_out_YM.Rdata")






